package scu.maqiang.fes;

import scu.maqiang.mesh.Mesh1L3;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

/**
 * 一维二次元有限元空间
 * @author 马强
 *
 */
public class FES1L31 extends FES{

	public FES1L31(Mesh1L3 mesh1D) {
		super(mesh1D, 1);
	}
	
	
	public void assembleStiff(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES1L31::elementStiff, tp, A);
	}
	
	public void assembleStiff(double[] constCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(constCoef, FES1L31::elementStiff, tp, A);
	}

	public void assembleStiff(double[][] vecCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(vecCoef, FES1L31::elementStiff, tp, A);
	}
	
	public static void elementStiff(int[] ele, double[][] xx, double[][] kxx, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xx, 0);
		double L = Math.abs(x[0] - x[1]);
		double coef;
		if (tp == BVPType.COMMON) {
			coef = (kxx[0][0] + kxx[0][1] + 4 * kxx[0][2]) / (6 * L);
		} else if (tp == BVPType.AXISYMMETRIC){
			coef = (kxx[0][0] * x[0] + kxx[0][1] * x[1] + 4 * kxx[0][2] * x[2]) / (6 * L);
		} else {
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		eleMatrix[0][0] = eleMatrix[1][1] = coef * 7.0/3;
		eleMatrix[0][1] = eleMatrix[1][0] = coef / 3;
		eleMatrix[0][2] = eleMatrix[2][0] = eleMatrix[1][2] = eleMatrix[2][1] = -coef * 8.0/3;
		eleMatrix[2][2] = coef * 16.0 / 3;
	}
	
	
	public void assembleHeatSource(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS) {
		assembleGlobalVector(coefFunc, param, FES1L31::elementSource, tp, RHS);
	}

	public void assembleHeatSource(double[] constCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(constCoef, FES1L31::elementSource, tp, RHS);
	}

	public void assembleHeatSource(double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(vecCoef, FES1L31::elementSource, tp, RHS);
	}

	public static void elementSource(int[] ele, double[][] xx, double[][] vf, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xx, 0);
		double L = Math.abs(x[1] - x[0]);
		double[] f = vf[0];
		double mf = 0.0;
		if (tp == BVPType.COMMON) {
			mf = (f[0] + f[1] + 4 * f[2]) / 6.0;
		} else if (tp == BVPType.AXISYMMETRIC) {
			mf = (f[0] * x[0] + f[1] * x[1] + 4 * f[2] * x[2]) / 6.0;
		} else if (tp == BVPType.SPHERICAL) {
			mf = (f[0] * x[0] * x[0] + f[1] * x[1] * x[1] + 4 * f[2] * x[2] * x[2]) / 6.0;
		}
		eleF[0] = eleF[1] = mf * L / 6;
		eleF[2] = 4 * mf * L / 6;
	}
	
	public void computeDerivatives(double[] u, double[] du, double[] ddu) {
	    double[] count = new double[ndof];
	    double[][] coord = new double[3][1];
	    double[] eleU = new double[3];
	    double h;
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	double[] x = MVO.col(coord, 0);
	    	eleU[0] = u[ele[0]];
	    	eleU[1] = u[ele[1]];
	    	eleU[2] = u[ele[2]];
	    	h = x[1] - x[0];
	    	du[ele[0]] += eleU[0] * (-3.0 / h) + eleU[1] * (-1.0 / h) + eleU[2] * ( 4.0 / h);
	    	du[ele[1]] += eleU[0] * ( 1.0 / h) + eleU[1] * ( 3.0 / h) + eleU[2] * (-4.0 / h);
	    	du[ele[2]] += (eleU[1] - eleU[0]) / h;
	    	ddu[ele[0]] = ddu[ele[1]] = ddu[ele[2]] = 4.0 / (h * h) *(eleU[0] - 2 * eleU[2] + eleU[1]);
	    	count[ele[0]] += 1.0;
	    	count[ele[1]] += 1.0;
	    	count[ele[2]] += 1.0;
	    }
	    for(int i = 0; i < ndof; i++) {
	    	du[i] /= count[i];
	    }
	}

}
